
# Minimal packages
import os, sys
from pystac_client import Client
from odc.stac import configure_s3_access, stac_load

# Python packages
import re
import json
import pandas as pd
import numpy as np
from pathlib import Path

# ODC packages
from dea_tools.plotting import display_map, rgb
from datacube.utils import masking
from odc.algo import mask_cleanup, erase_bad, to_f32
from odc.ui import image_aspect

# EASI packages
repo = Path.home() / 'eocsi-hackathon-2022'  # No easy way to get repo directory
if repo not in sys.path: sys.path.append(str(repo))
from tools.notebook_utils import heading, xarray_object_size, initialize_dask, localcluster_dashboard
from tools.stac_utils import stac_landsat_assets_df, stac_landsat_flags_to_dc


# Setup

# Does this work stand-alone or require an AWS account?
configure_s3_access(requester_pays=True)

# Optional: use EASI SE Asia caching-proxy service
os.environ["AWS_HTTPS"] = "NO"
os.environ["GDAL_HTTP_PROXY"] = "easi-caching-proxy.caching-proxy:80"
print(f'Will use caching proxy at: {os.environ.get("GDAL_HTTP_PROXY")}')

# Optional: Dask
cluster, client = initialize_dask(use_gateway=True, workers=(1,7), wait=False)
if cluster: display(cluster)
if cluster is None or 'LocalCluster' in str(type(cluster)): display(localcluster_dashboard(client))

Will use caching proxy at: easi-caching-proxy.caching-proxy:80
An existing cluster was found. Connecting to: easihub.1df2f465d26344ed9e723be1b41f5d75

VBox(children=(HTML(value='<h2>GatewayCluster</h2>'), HBox(children=(HTML(value='\n<div>\n<style scoped>\n    …


# Vietnam - Ha Long
latitude = (20.5, 21.1)
longitude = (106.5, 107.2)
time=('2022-01-01', '2022-06-01')

# PNG Milne Bay
# latitude = (-10.8, -10)
# longitude = (149.7, 150.8)  
# time=('2022-01-01', '2022-03-01')

# Fiji - blows up JHub memory due to antemeridian
# latitude = (-17.1, -16.2)
# longitude = (178.2, 180.0)
# time=('2020-02-01', '2020-02-20')

# west, south, east, north
bbox = [longitude[0], latitude[0], longitude[1], latitude[1]]

# Display bounding box on a map
display_map(longitude, latitude)


if do_landsat:

    # STAC catalog and query
    catalog = Client.open('https://landsatlook.usgs.gov/stac-server/')
    product = 'landsat-c2l2-sr'
    query_cfg = ["platform=LANDSAT_8", "landsat:collection_category=T1"]

    # Search for available items
    query = catalog.search(
        collections=[product], datetime=f'{time[0]}/{time[1]}', bbox=bbox, query=query_cfg
    )
    items = list(query.get_items())
    print(f"Found: {len(items):d} datasets")

    # Rewrite URLs to use S3
    def landsat_patch(uri: str) -> str:
        """Return the S3 version of the URI"""
        return uri.replace('https://landsatlook.usgs.gov/data/', 's3://usgs-landsat/')

    # Change or update STAC information for use by ODC
    stac_cfg = json.load(open(f'{repo}/stac_cfgs/aws/{product}.json'))

    # `stac_load` parameters
    bands = ('blue', 'green', 'red', 'nir08', 'qa_pixel')
    stac_call = {
        'bands': bands,                            # Optional: selected bands
        'bbox': bbox,                              # Bounding box. Also Geopolygon, GeoBox or None (full extent of items)
        'chunks': {'x': 4096, 'y': 4096},          # Optional: if using Dask
        'groupby': "solar_day",                    # "solar_day" = group scenes on same solar day into same time layer in cube
        'stac_cfg': stac_cfg,
        'patch_url': landsat_patch,
    }

    # Additional Landsat band specifications
    band_specs = {
        'red': {
            'scale': 0.0000275,
            'offset': -0.2
        },
        'nir08': {
            'scale': 0.0000275,
            'offset': -0.2
        },
    }

Found: 10 datasets


if do_sentinel:
    
    # STAC catalog and query
    catalog = Client.open('https://earth-search.aws.element84.com/v0')
    product = 'sentinel-s2-l2a-cogs'
    
    # Search for available items
    query = catalog.search(
        collections=[product], datetime=f'{time[0]}/{time[1]}', bbox=bbox,
    )
    items = list(query.get_items())
    print(f"Found: {len(items):d} datasets")
    
    # Rewrite URLs to use S3
    def patch(uri: str) -> str:
        """Return the Sentinel-2 S3 version of the URI"""
        return uri.replace('https://sentinel-cogs.s3.us-west-2.amazonaws.com/', 's3://sentinel-cogs/')
    
    # Change or update STAC information for use by ODC 
    stac2odc_cfg = {
        "sentinel-s2-l2a-cogs": {
            "assets": {
                "*": {"data_type": "uint16", "nodata": 0},
                "SCL": {"data_type": "uint8", "nodata": 0},
                "visual": {"data_type": "uint8", "nodata": 0},
            },
            "aliases": {"red": "B04", "green": "B03", "blue": "B02"},
        },
        "*": {"warnings": "ignore"},
    }
    
    # `stac_load` parameters
    stac_call = {
        'bands': ("B04",),
        'crs': crs,
        'resolution': 30,
        # chunks={},  # <-- use Dask
        # groupby="solar_day",
        'stac_cfg': cfg,
        'patch_url': patch,
    }


# Optional: Explore the structure of a STAC item
display(items[0])

# Optional: List available band names and selected details
display(stac_landsat_assets_df(items[0]))


xx = stac_load(items, **stac_call)

heading(xarray_object_size(xx))
display(xx)
display(xx.odc.geobox)
aspect = image_aspect(xx)

<xarray.Dataset>
Dimensions:      (y: 2243, x: 2455, time: 6)
Coordinates:
  * y            (y) float64 2.335e+06 2.335e+06 ... 2.268e+06 2.268e+06
  * x            (x) float64 6.558e+05 6.558e+05 ... 7.294e+05 7.294e+05
    spatial_ref  int32 32648
  * time         (time) datetime64[ns] 2022-01-02T03:17:29.057704 ... 2022-05...
Data variables:
    blue         (time, y, x) uint16 dask.array<chunksize=(1, 2243, 2455), meta=np.ndarray>
    green        (time, y, x) uint16 dask.array<chunksize=(1, 2243, 2455), meta=np.ndarray>
    red          (time, y, x) uint16 dask.array<chunksize=(1, 2243, 2455), meta=np.ndarray>
    nir08        (time, y, x) uint16 dask.array<chunksize=(1, 2243, 2455), meta=np.ndarray>
    qa_pixel     (time, y, x) uint16 dask.array<chunksize=(1, 2243, 2455), meta=np.ndarray>

array([2334795., 2334765., 2334735., ..., 2267595., 2267565., 2267535.])

array([655815., 655845., 655875., ..., 729375., 729405., 729435.])

array(32648, dtype=int32)

array(['2022-01-02T03:17:29.057704000', '2022-01-18T03:17:26.701116000',
       '2022-03-07T03:17:14.069855000', '2022-04-08T03:17:01.924989000',
       '2022-05-10T03:17:31.590688000', '2022-05-26T03:17:09.908311000'],
      dtype='datetime64[ns]')


# Convert pixel quality flag descriptions to ODC format
flags_def = stac_landsat_flags_to_dc(
    items[0].assets['qa_pixel'].to_dict().get('classification:bitfields')
)
# heading('"qa_pixel" flags definition')
# display(flags_def)

# Cloud mask flags (logical AND?)
quality_flags = {
    'cloud': 'cloud',    # True where there is cloud
    'cirrus': 'cirrus',  # True where there is cirrus cloud
    'shadow': 'shadow',  # True where there is cloud shadow 
}

# Set bit mask: True=cloud, False=non-cloud
mask, _= masking.create_mask_value(flags_def, **quality_flags)

# Add the cloud mask to our dataset
xx['cloud_mask'] = (xx['qa_pixel'] & mask) != 0  # bitwise-and != 0 simulates an 'or'

# Apply morphological processing on the cloud mask
filters = [("opening", 2),("dilation", 2)]
xx['cloud_mask_filtered'] = mask_cleanup(xx['cloud_mask'], mask_filters=filters)

# Apply the the cloud-mask to the data variables
clear_filtered = erase_bad(xx.drop_vars(['cloud_mask_filtered', 'cloud_mask', 'qa_pixel']),
                           xx['cloud_mask_filtered'])

# Apply scale and offset
scale = stac_cfg[product]['assets']['*']['scale']
offset = stac_cfg[product]['assets']['*']['offset']
clear_filtered = (clear_filtered * scale + offset).astype(np.float32)


xx = xx.persist()

clear_filtered = clear_filtered.persist()


# Optional: Check the data structures and sizes

# heading(xarray_object_size(xx))
# display(xx)

# heading(xarray_object_size(clear_filtered))
# display(clear_filtered)


# Optional: Plot and RGB timeseries
rgb(xx, ['red', 'green', 'blue'], col='time', col_wrap=4)

# Optional: plot the masked and scaled data
rgb(clear_filtered, ['red', 'green', 'blue'], col='time', col_wrap=4)

2022-09-29 08:52:37,154 - distributed.client - ERROR - Failed to reconnect to scheduler after 30.00 seconds, closing client


# client.close()
# if cluster: cluster.close()

	title	description	eo:bands	classification:bitfields
coastal	Coastal/Aerosol Band (B1)	Collection 2 Level-2 Coastal/Aerosol Band (B1)...	[{'name': 'B1', 'common_name': 'coastal', 'gsd...	NaN
blue	Blue Band (B2)	Collection 2 Level-2 Blue Band (B2) Surface Re...	[{'name': 'B2', 'common_name': 'blue', 'gsd': ...	NaN
green	Green Band (B3)	Collection 2 Level-2 Green Band (B3) Surface R...	[{'name': 'B3', 'common_name': 'green', 'gsd':...	NaN
red	Red Band (B4)	Collection 2 Level-2 Red Band (B4) Surface Ref...	[{'name': 'B4', 'common_name': 'red', 'gsd': 3...	NaN
nir08	Near Infrared Band 0.8 (B5)	Collection 2 Level-2 Near Infrared Band 0.8 (B...	[{'name': 'B5', 'common_name': 'nir08', 'gsd':...	NaN
swir16	Short-wave Infrared Band 1.6 (B6)	Collection 2 Level-2 Short-wave Infrared Band ...	[{'name': 'B6', 'common_name': 'swir16', 'gsd'...	NaN
swir22	Short-wave Infrared Band 2.2 (B7)	Collection 2 Level-2 Short-wave Infrared Band ...	[{'name': 'B7', 'common_name': 'swir22', 'gsd'...	NaN
qa_aerosol	Aerosol Quality Analysis Band	Collection 2 Level-2 Aerosol Quality Analysis ...	NaN	[{'name': 'fill', 'description': 'Correspondin...
qa_pixel	Pixel Quality Assessment Band	Collection 2 Level-2 Pixel Quality Assessment ...	NaN	[{'name': 'fill', 'description': 'Correspondin...
qa_radsat	Radiometric Saturation Quality Assessment Band	Collection 2 Level-2 Radiometric Saturation Qu...	NaN	[{'name': 'band1', 'description': 'Band 1 radi...

Download Landsat data using STAC API¶

Landsat data¶

Data source and documentation¶

Spatio Temporal Asset Catalogs (STAC)¶

Open Data Cube¶

More information¶

Select an area of interest¶

Search USGS STAC catalog¶

Sentinel-2 configuration and settings¶

Item: LC08_L2SP_126046_20220526_20220602_02_T1_SR

STAC Extensions

Assets

Asset: Thumbnail image

Asset: Reduced resolution browse image

Asset: HTML index page

Asset: Product Metadata File (json)

Asset: Coastal/Aerosol Band (B1)

Asset: Blue Band (B2)

Asset: Green Band (B3)

Asset: Red Band (B4)

Asset: Near Infrared Band 0.8 (B5)

Asset: Short-wave Infrared Band 1.6 (B6)

Asset: Short-wave Infrared Band 2.2 (B7)

Asset: Aerosol Quality Analysis Band

Asset: Pixel Quality Assessment Band

Asset: Radiometric Saturation Quality Assessment Band

Asset: Angle Coefficients File

Asset: Product Metadata File

Asset: Product Metadata File (xml)

Links

Link:

Link:

Link:

Link:

Load the selected items into an xarray cube¶

Dataset size: 315.12 MB

GeoBox

Apply pixel quality masking and scaling¶

Dask: compute results¶

Optional plot summaries¶

Export to netCDF or Geotiff¶

Close the Dask cluster¶

Load the selected items into an `xarray` cube¶